//
// CUDA Sagan Simulator - a CUDA version fo the sagan Simulator Software
// version 1.0 - Alpha 01
// by TTYMan
//
// para compilar:  nvcc SS_v3.cu -o sagan -I/home/ttyman/NVIDIA_GPU_Computing_SDK/C/common/inc -arch sm_21

#include <stdio.h>
#include <string.h>
#include <cutil.h>
#include <cuda.h>
#include <cuda_runtime.h>

/*
#if __CUDA_ARCH__ < 130
typedef float REAL;
#else
typedef double REAL;
#endif
*/

typedef float  REAL;
typedef float4 REAL_CUDA;

//cutil.h ==> /home/ttyman/NVIDIA_GPU_Computing_SDK/C/common/inc

#define BSIZE 256
#define BUFFER_MAX 2048      // buffer size for reading each line of initial conditions file.
#define offset      4        // this offset space is used to define particle position information (float4 CUDA datatype):
                             //   0 (x) - Pos X ==> array[index*offset+0]
                             //   1 (y) - Pos Y ==> array[index*offset+1]
                             //   2 (z) - Pos Z
                             //   3 (w) - Mass
                             // The velocity information is stored ina float3 CUDA datatype, as follows:
                             //   0 (x) - Vel X
                             //   1 (y) - Vel Y
                             //   2 (z) - Vel Z
                             //   3 (W) - Body color (to be used in future implementations, not used now).
                             // this approuch is useful to take advantage of CUDA float3 and float4 datatypes.


#define SHARED_KERNEL 96
#define SHARED_MEMORY_SIZE 4096

#if __CUDA_ARCH__ <300
#define BLOCKS_PER_MP 8
#else
#define BLOCKS_PER_MP 16
#endif

// declaracoes memoria device

// device memory for particles array
__device__ REAL_CUDA*  d_ParticleArrayPosition;
__device__ REAL_CUDA*  d_ParticleArrayVelocity;
int quantidade;


// Kernels
__device__ REAL_CUDA BodyBodyInteraction(REAL_CUDA ParticleIPosition, REAL_CUDA ParticleJPosition, REAL_CUDA Acceleration, REAL SofteningFactorSQR)
{
    REAL XDist;
    REAL YDist;
    REAL ZDist;
    REAL DistanceSQR;
    REAL Rinv;
    REAL Rinv3;

	// [1 FLOP]
	XDist = ParticleIPosition.x - ParticleJPosition.x;

	// [1 FLOP]
    YDist = ParticleIPosition.y - ParticleJPosition.y;

    // [1 FLOP]
	ZDist = ParticleIPosition.z - ParticleJPosition.z;

    // [5 FLOP]
	DistanceSQR = (XDist * XDist) + (YDist * YDist) + (ZDist * ZDist);

    // [4 FLOP]
    Rinv = 1.0 / sqrt( DistanceSQR + SofteningFactorSQR ); 

    // [2 FLOP]
    Rinv3 = (Rinv * Rinv * Rinv);

    // [9 FLOP]
    Acceleration.x -= ParticleJPosition.w * XDist * Rinv3;
	Acceleration.y -= ParticleJPosition.w * YDist * Rinv3;
	Acceleration.z -= ParticleJPosition.w * ZDist * Rinv3;

    return Acceleration;
	

}

__device__ void FillSharedMemoryBuffer(REAL_CUDA* d_ParticleArrayPosition, unsigned int N, unsigned long tile, unsigned long tilesAmount, int tileWidth)
{
	extern __shared__ REAL_CUDA sh_ParticleArrayPosition[];
	
	// last tile?
	if(tile == tilesAmount-1)
	{
		for(unsigned int i=0; i<tileWidth; i++)
		{
			if(tile+i < N)
			{
				sh_ParticleArrayPosition[i] = d_ParticleArrayPosition[tile*tileWidth+i];
				/*sh_ParticleArrayPosition[i].x = 1.0;
				sh_ParticleArrayPosition[i].y = 1.0;
				sh_ParticleArrayPosition[i].z = 1.0;
				sh_ParticleArrayPosition[i].w = 1.0;*/
			} else break; // avoid out of bound access on shared memory buffer.
		}		
	}
	else
	{
		for(unsigned int i=0; i<tileWidth; i++)
		{
			sh_ParticleArrayPosition[i] = d_ParticleArrayPosition[tile*tileWidth+i];
			/*sh_ParticleArrayPosition[i].x = 1.0;
			sh_ParticleArrayPosition[i].y = 1.0;
			sh_ParticleArrayPosition[i].z = 1.0;
			sh_ParticleArrayPosition[i].w = 1.0;*/
		}
	}
	
}

__device__ REAL_CUDA TileProcessing(REAL_CUDA* d_ParticleArrayPosition, unsigned int N, REAL_CUDA Particle, REAL_CUDA Acceleration, REAL SofteningFactorSQR, unsigned long tile, unsigned long tilesAmount, int tileWidth)
{
	extern __shared__ REAL_CUDA sh_ParticleArrayPosition[];
	
	// last tile?
	if(tile == tilesAmount-1)
	{
		for(unsigned int i=0; i<tileWidth; i++)
		{
			if(tile+i < N)
			{
				Acceleration = BodyBodyInteraction(Particle, sh_ParticleArrayPosition[i], Acceleration, SofteningFactorSQR);
			} else break; // avoid out of bound access on shared memory buffer.
		}		
	}
	else
	{
		for(unsigned int i=0; i<tileWidth; i++)
		{
			Acceleration = BodyBodyInteraction(Particle, sh_ParticleArrayPosition[i], Acceleration, SofteningFactorSQR);
		}
	}
	
	// teste teste
	Acceleration.x = 0.0;
	Acceleration.y = 0.0;
	Acceleration.z = 0.0;
	Acceleration.w = 0.0;
	return Acceleration;
}


__global__ void MainKernel(REAL_CUDA* d_ParticleArrayPosition,
                           REAL_CUDA* d_ParticleArrayVelocity,
                           unsigned int N,
                           REAL SofteningFactor,
                           REAL TimeStep,
                           REAL	*d_flag, // para que serve isso, Lucas??? Vou manter por enquanto.
                           int tileWidth)
{
	// map shared memory buffer
	extern __shared__ REAL_CUDA sh_ParticleArrayPosition[];
	
	// global thread ID: unique ID for the particle.
	unsigned int gThID = blockIdx.x * blockDim.x + threadIdx.x; 
	
	REAL_CUDA Particle = d_ParticleArrayPosition[gThID]; //global memory access!
	
	REAL_CUDA Acceleration = {0.0, 0.0, 0.0, 0.0};
	
	// save some FLOPs squaring epson factor here, avoiding to do this inside the tile processing loop.
	REAL SofteningFactorSQR = SofteningFactor * SofteningFactor;
		
	// how many tiles will be used?
	unsigned long tilesAmount = N / tileWidth;
	//int tilesAmount = 12;

	for(unsigned long tile=0; tile < tilesAmount; tile++)
	{

		// the first block thread will be used to fill the shared memory buffer.
		if(threadIdx.x==0)
		{
			FillSharedMemoryBuffer(d_ParticleArrayPosition, N, tile, tilesAmount, tileWidth);
			
		}
		__syncthreads(); //make sure shared memory buffer was updated for all threads.
		
		Acceleration = TileProcessing(d_ParticleArrayPosition, N, Particle, Acceleration, SofteningFactorSQR, tile, tilesAmount, tileWidth);

	}
	
	
	//after processing all tiles, saves back the acceleration value on the global memory array.
	//REAL_CUDA accTemp = {Acceleration.x, Acceleration.y, Acceleration.z, Acceleration.w};
	//d_ParticleArrayVelocity[gThID] = Acceleration; //accTemp; // global memory access!
	
	__syncthreads();
	
	// updating the velocities
	d_ParticleArrayVelocity[gThID].x   += Acceleration.x * TimeStep;
	d_ParticleArrayVelocity[gThID].y   += Acceleration.y * TimeStep;
	d_ParticleArrayVelocity[gThID].z   += Acceleration.z * TimeStep;

	// updating the coordinates
	d_ParticleArrayPosition[gThID].x += (d_ParticleArrayVelocity[gThID].x * TimeStep)/2.0;
	d_ParticleArrayPosition[gThID].y += (d_ParticleArrayVelocity[gThID].y * TimeStep)/2.0;
	d_ParticleArrayPosition[gThID].z += (d_ParticleArrayVelocity[gThID].z * TimeStep)/2.0;	
	
}





void EnergyLog(char* FileName, REAL time, REAL Etot, REAL Ek, REAL Ep, REAL energyerr, char* fmode)
{

    //
    // used to generated the Energy Log.
    //

    FILE *file_p;

    file_p = fopen(FileName, fmode);

	fprintf(file_p, "%e %e %e %e %e\n", time, Etot, Ek, Ep, energyerr);

    fclose(file_p);

}


void LoadInitialConditionsData(char* FileName, unsigned int* N, REAL_CUDA** ParticleArrayPosition, REAL_CUDA** ParticleArrayVelocity)
{
    FILE* fptr;

    if ((fptr = fopen(FileName, "r")))
    {
        char    LineBuffer[BUFFER_MAX];
        unsigned int     LineCount = 0;

        // The first line of the initial conditions file should contain the amount of particles information.
        // So, read this value and use it for allocate memory for the particles arrays.
        fgets (LineBuffer, BUFFER_MAX, fptr);
        sscanf(LineBuffer, "%d ", &LineCount);
        *N = LineCount; // 'N' value finally discovered!

        // allocate memory for the Particle arrays.
        *ParticleArrayPosition = (REAL_CUDA*)malloc( sizeof(REAL_CUDA) * LineCount );
        *ParticleArrayVelocity = (REAL_CUDA*)malloc( sizeof(REAL_CUDA) * LineCount );

        // read the file line by line, saving position and velocity into data arrays.
        // each file line represents only one particle.
        for (int i=0; i<LineCount; i++)
        {
            fgets (LineBuffer, BUFFER_MAX, fptr);
            //sscanf(LineBuffer, "%lf %lf %lf %lf %lf %lf %lf " , &(*ParticleArrayPosition)[i].x,      // PosX ==> &(((*ParticleArrayPosition)+i)->x),
            sscanf(LineBuffer, "%f %f %f %f %f %f %f " , &(*ParticleArrayPosition)[i].x,      // PosX ==> &(((*ParticleArrayPosition)+i)->x),
                                                         &(*ParticleArrayPosition)[i].y,      // PosY
                                                         &(*ParticleArrayPosition)[i].z,      // PosZ
                                                         &(*ParticleArrayVelocity)[i].x,      // VelX
                                                         &(*ParticleArrayVelocity)[i].y,      // VelY
                                                         &(*ParticleArrayVelocity)[i].z,      // VelZ
                                                         &(*ParticleArrayPosition)[i].w);     // Mass
        }

        //close the file and release ay used resource
        fclose(fptr);

        return;

    }
    else
    {
        //Error! Couldn't open the file!
        printf("\nError opening file %s!\n", FileName);
        return;
    }
}

void GenerateResultFile(char* FileName, REAL_CUDA* ParticleArrayPosition, REAL_CUDA* ParticleArrayVelocity, unsigned int N)
{

    //
    // used to generated the data file containing the updated particles positions.
    //

    FILE *file_p;

    file_p = fopen(FileName, "w");

    //first line must cointain the number of particles
    //fprintf(file_p, "%d\n", N);

    for(int i=0; i<N; i++)
    {
        //fprintf(file_p, "%f %f %f %f %f %f %f\n", ParticleArrayPosition[i].x, // teste com uso de double4
	    fprintf(file_p, "%e %e %e %e %e %e %e\n", ParticleArrayPosition[i].x,
                                                  ParticleArrayPosition[i].y,
                                                  ParticleArrayPosition[i].z,
                                                  ParticleArrayVelocity[i].x,
                                                  ParticleArrayVelocity[i].y,
                                                  ParticleArrayVelocity[i].z,
                                                  ParticleArrayPosition[i].w);

    }

    fclose(file_p);

}

void traceParticle(int id, REAL t, REAL* KE, REAL* PE, char* filemode)
{
	char ResultFileName[100];

	sprintf(ResultFileName, "./data/energy-particule%d.log", id);

	FILE *file_p;

    file_p = fopen(ResultFileName, filemode);

	fprintf(file_p, "%e %e %e\n", t, *KE, *PE);

    fclose(file_p);

    //soma das energias
    sprintf(ResultFileName, "./data/somaEnergy-particule%d.log", id);
    file_p = fopen(ResultFileName, filemode);

	fprintf(file_p, "%e\n",(REAL)(*KE+*PE));

    fclose(file_p);
}

int main(int argc, char **argv)
{
    printf("\nCUDASaganSimulator version 1.0 Beta\nStarting...\n:: Getting information about available GPU...\n");


	// Getting information about installed GPU
	// this information will be used to calculate the size of execution kernel blocks.
    int devID;
    int MAX_NUMBER_THREADS_PER_BLOCK;
    int NUMBER_OF_BLOCKS;
    int THREADS_PER_BLOCK;
    cudaDeviceProp props;
    CUDA_SAFE_CALL( cudaGetDevice(&devID) );
    CUDA_SAFE_CALL( cudaGetDeviceProperties(&props, devID) );
    if (props.major == 9999 && props.minor == 9999)
    {
        // there is not CUDA devices available!!!
        printf("\n ==>>> There is no CUDA device available! Exiting now...\n");
        return 0;
    }
    MAX_NUMBER_THREADS_PER_BLOCK = props.maxThreadsPerBlock;

    printf(":::: GPU found: %s\n", props.name);
    printf(":::: Maximum Number os threads per block allowed: %d\n", MAX_NUMBER_THREADS_PER_BLOCK);




    // discrete variables
    //char InitialConditionsFileName[50];
    char InitialConditionsFileName[] = "esfera395CUDA.ini";
    char strTrace[20];
    int trace;
    char    Arquivo[100];
    unsigned int     N = 0;

    if (argc>=2)
    {
		int argcont;

		for(argcont=0;argcont<argc;argcont++)
		{

			if (!strcmp(argv[argcont],"-init"))
			{
				strcpy(InitialConditionsFileName,argv[argcont+1]);
			}

			if (!strcmp(argv[argcont],"-trace"))
			{
				strcpy(strTrace,argv[argcont+1]);
			}
		}
	}

	else
	{
		printf("\n:::Initial file not given, please start the program using -init filename.ini:::\n\n");
		exit(1);
	}

	sprintf(Arquivo,"./data/%s",InitialConditionsFileName);
	if (strlen(strTrace)>0)
		trace = atoi(strTrace);
	else
		trace = 0;

	// trace the "id" of particle
	REAL* traceKE;
	REAL* tracePE;
	REAL* d_traceKE;
	REAL* d_tracePE;
	traceKE   = (REAL*)malloc(sizeof(REAL));
    tracePE = (REAL*)malloc(sizeof(REAL));
    cudaMalloc( (void**)&d_traceKE,   sizeof(REAL));
    cudaMalloc( (void**)&d_tracePE, sizeof(REAL));

    REAL flag = 0;
    REAL* d_flag;


    // host memory for particle arrays.
    REAL_CUDA* ParticleArrayPosition;
    REAL_CUDA* ParticleArrayVelocity;


    // Load 'N' value and particle information from initial conditions file.
    printf(":: Loading Initial Conditions File (%s)...\n", Arquivo);
    LoadInitialConditionsData(Arquivo, &N, &ParticleArrayPosition, &ParticleArrayVelocity);

    // Once 'N' is now knowed, it's possible to define the size of the execution block
    printf(":: Initial Conditions File loaded. %d particles were found.\n", N);

    /*
    if(N>MAX_NUMBER_THREADS_PER_BLOCK)
    {
        NUMBER_OF_BLOCKS =  ((N % MAX_NUMBER_THREADS_PER_BLOCK)>0 ? 1 : 0) + (N / MAX_NUMBER_THREADS_PER_BLOCK);
        THREADS_PER_BLOCK = MAX_NUMBER_THREADS_PER_BLOCK;
    }
    else
    {
        NUMBER_OF_BLOCKS = 1;
        THREADS_PER_BLOCK = N;
    }
    */
	float cc;
	cc = props.major+props.minor*0.1;
   	printf("\nCC: %f",cc);

	int THREADS_KERNEL;

	if (cc <1.3)
		THREADS_KERNEL = 16;
	else if (cc==1.3)
			THREADS_KERNEL = 26;
	else
			THREADS_KERNEL = 34;

	printf("\nThreads kernel: %d",THREADS_KERNEL);

	int warps_block = floor(props.regsPerBlock/(THREADS_KERNEL*props.warpSize));
	THREADS_PER_BLOCK = warps_block * props.warpSize;

	if (THREADS_PER_BLOCK>MAX_NUMBER_THREADS_PER_BLOCK)
	{
		THREADS_PER_BLOCK = MAX_NUMBER_THREADS_PER_BLOCK;
	}

	//este vai ser o numero maximo sem estourar o numero de regs
	MAX_NUMBER_THREADS_PER_BLOCK = THREADS_PER_BLOCK;

	int it;
	int warps_per_sm;
	int warps_per_block;
	int regs_per_block;
	int blocks_per_sm_regs;
	int blocks_per_sm_warps;
	int blocks_per_sm_shared;
	int warps_alloc_granularity;
	int rounded_warps_per_block;
	int regs_alloc_unit_size;
	int max_warps_per_sm;
	int max_blocks_per_sm;

	//set warps_alloc_granularity, values from nvidia
	if (cc >= 3.0)
	{
		warps_alloc_granularity = 4;
		max_blocks_per_sm = 16;
	}
	else
	{
		warps_alloc_granularity = 2;
		max_blocks_per_sm = 8;
	}

	//set regs_alloc_unit_size and max_warps_per_sm, values from nvidia
	if (cc < 1.2)
	{
		regs_alloc_unit_size = 256;
		max_warps_per_sm = 24;
	}
	else if (cc < 2)
	{
		regs_alloc_unit_size = 512;
		max_warps_per_sm = 32;
	}
	else if (cc< 3)
	{
		regs_alloc_unit_size = 64;
		max_warps_per_sm = 48;
	}
	else
	{
		regs_alloc_unit_size = 256;
		max_warps_per_sm = 64;
	}



	//max_threads é o numero de threads por bloco com maior ocupação
	//max_warps_per_mp é o numero de warps que o max_threads vai ocupar do mp
	int best_warps_per_sm=0;
	int best_threads=0;

	/*ARRUMAR REGS_PER_BLOCK */
	for(it=32;it<MAX_NUMBER_THREADS_PER_BLOCK;it+=32)
	{
		warps_per_block = it/props.warpSize;

		if (cc<2.0) //so allocation granularity = block
		{
			printf("\nallocation granularity = block");
			if (warps_per_block%warps_alloc_granularity != 0) //arredondar para cima warps_per_block
			{
				rounded_warps_per_block=(warps_per_block/warps_alloc_granularity)*warps_alloc_granularity+warps_alloc_granularity;
				regs_per_block=rounded_warps_per_block*THREADS_KERNEL*props.warpSize;
				if (regs_per_block%regs_alloc_unit_size != 0) //arredondar para cima regs_per_block
					regs_per_block = (regs_per_block/regs_alloc_unit_size)*regs_alloc_unit_size+regs_alloc_unit_size;
			}
		}
		else //so allocation granularity = warps
		{
			printf("\nallocation granularity = warps");
			if (THREADS_KERNEL*props.warpSize%regs_alloc_unit_size != 0) //arredondar para cima THREADS_KERNEL*props.warpSize
			{
				if (warps_per_block%warps_alloc_granularity != 0) //arredondar para cima warps_per_block
					regs_per_block = ((THREADS_KERNEL*props.warpSize/regs_alloc_unit_size)*regs_alloc_unit_size+regs_alloc_unit_size)*((warps_per_block/warps_alloc_granularity)*warps_alloc_granularity+warps_alloc_granularity);
				else
					regs_per_block = ((THREADS_KERNEL*props.warpSize/regs_alloc_unit_size)*regs_alloc_unit_size+regs_alloc_unit_size)*warps_per_block;
			}
			else //nao precisei arrendondar THREADS_KERNEL*props.warpSize
			{
				//OK
				if (warps_per_block%warps_alloc_granularity != 0) //arredondar para cima warps_per_block
					regs_per_block = THREADS_KERNEL*props.warpSize*((warps_per_block/warps_alloc_granularity)*warps_alloc_granularity+warps_alloc_granularity);
				else //nao precisei arrendondar warps_per_block
					regs_per_block = THREADS_KERNEL*props.warpSize*warps_per_block;
			}
		}

		//blocos por sm devido ao numero de regs
		blocks_per_sm_regs = props.regsPerBlock/regs_per_block;
		//printf("\nprops.regsPerBlock: %d, regs_per_block: %d, blocks_per_sm_regs: %d",props.regsPerBlock,regs_per_block,blocks_per_sm_regs);

		//blocos por sm devido ao numero de warps
		blocks_per_sm_warps = max_warps_per_sm/warps_per_block;
		if (blocks_per_sm_warps>max_blocks_per_sm) //minimo deles
			blocks_per_sm_warps=max_blocks_per_sm;

		//blocos por sm devido ao numero de shared memory
		blocks_per_sm_shared = props.sharedMemPerBlock/SHARED_MEMORY_SIZE;

		if (blocks_per_sm_regs <= blocks_per_sm_warps && blocks_per_sm_regs <= blocks_per_sm_shared)
			warps_per_sm = blocks_per_sm_regs * warps_per_block;
		else if (blocks_per_sm_warps <= blocks_per_sm_shared)
			warps_per_sm = blocks_per_sm_warps * warps_per_block;
		else
			warps_per_sm = blocks_per_sm_shared * warps_per_block;

		//printf("\nblocks_per_sm_regs: %d, blocks_per_sm_warps: %d, blocks_per_sm_shared: %d, warps_per_sm: %d",blocks_per_sm_regs,blocks_per_sm_warps,blocks_per_sm_shared,warps_per_sm);

		if (warps_per_sm > best_warps_per_sm && (N/it + N%it>0?1:0)<=props.multiProcessorCount*max_blocks_per_sm)
		{
			best_warps_per_sm = warps_per_sm;
			best_threads = it;
		}
	}

	THREADS_PER_BLOCK = best_threads;

	//teste
	//THREADS_PER_BLOCK = 320;

	if (N <= THREADS_PER_BLOCK)
	{
		NUMBER_OF_BLOCKS = 1;
	}

	else
	{
		NUMBER_OF_BLOCKS =  (N % THREADS_PER_BLOCK > 0 ? 1 : 0) + N / THREADS_PER_BLOCK;
	}


    printf(":: Number of execution blocks = %d\n", NUMBER_OF_BLOCKS);
    printf(":: Number of threads per block = %d\n", THREADS_PER_BLOCK);



    // device memory allocation for particle data arrays.
    CUDA_SAFE_CALL ( cudaMalloc((void**)&d_ParticleArrayPosition, sizeof(REAL_CUDA)*N) );
    if ( cudaSuccess != cudaGetLastError() )
	    printf( "Error malloc position!\n" );

    CUDA_SAFE_CALL ( cudaMalloc((void**)&d_ParticleArrayVelocity, sizeof(REAL_CUDA)*N) );
    if ( cudaSuccess != cudaGetLastError() )
	    printf( "Error malloc velocity!\n" );

    if (d_ParticleArrayPosition ==NULL || d_ParticleArrayVelocity == NULL)
    {
		printf("\nErro no cudaMalloc!");
		exit(1);
	}


    // copy particle information already loaded from initial conditions file to data arrays on device memory.
    CUDA_SAFE_CALL(cudaMemcpy(d_ParticleArrayPosition, ParticleArrayPosition, sizeof(REAL_CUDA)*N, cudaMemcpyHostToDevice));
    if ( cudaSuccess != cudaGetLastError() )
	{
		printf("erro no cudamemcpy");
	}

    CUDA_SAFE_CALL( cudaMemcpy(d_ParticleArrayVelocity, ParticleArrayVelocity, sizeof(REAL_CUDA)*N, cudaMemcpyHostToDevice) );
	if ( cudaSuccess != cudaGetLastError() )
	{
		printf("erro no cudamemcpy");
	}

    // Run the Kernel => no blocks, only threads for now :-)
    //UpdateParticleData<<<1, N>>>(d_ParticleArrayPosition, d_ParticleArrayVelocity);

    int i=0;    // we need these sequential numbers to name the output files.
    REAL t = 0.0;        // controls the amount of time steps already processed for the simulation.
    REAL dt = 0.0001;    // indicates the contribution of each scan for the simulation evolution.
    REAL dtout = 0.01;   // indicates the value of the step where a result file has to be generated.
    REAL tout = t+dtout;
    int    timeSteps = 10; //WARNING!! Don't forget to update this value with the desired time steps amount!
    REAL epson = 0.005;
    if(N>100)epson = 4.0/N; // softening factor

    REAL somErr = 0;
    int contErr = 0;

    char ResultFileName[100];

    // Energy calculation
    REAL* KineticEnergy;
    REAL* PotentialEnergy;
    REAL* d_KineticEnergy;
    REAL* d_PotentialEnergy;
    KineticEnergy   = (REAL*)malloc(sizeof(REAL));
    PotentialEnergy = (REAL*)malloc(sizeof(REAL));
    cudaMalloc( (void**)&d_KineticEnergy,   sizeof(REAL));
    cudaMalloc( (void**)&d_PotentialEnergy, sizeof(REAL));

    //flag
    CUDA_SAFE_CALL(cudaMalloc( (void**)&d_flag,sizeof(REAL)));
    CUDA_SAFE_CALL( cudaMemcpy(d_flag, &flag, sizeof(REAL), cudaMemcpyHostToDevice) );


	//
	// Definindo o tamanho da sharedmemory a ser alocada. queremos limitar em 16Kb por bloco
	//
	int tileWidth = 256; // quantidade de posiçoes a serem alocadas na shared memory por bloco. 
	// assumindo sizeof(double) = 16bytes, logo double4 = 16 * 4 = 64bytes
	// entao, 256 posicoes de double4 = 16K
	size_t SharedMemorySize = tileWidth * sizeof(REAL_CUDA);
	
    //printf(":: Number of execution blocks = %d\n", NUMBER_OF_BLOCKS);
    //printf(":: Number of threads per block = %d\n", THREADS_PER_BLOCK);	
	
	
	

    // this first execution is used to calculate the initial energy of the system. No updates are done.
    //size_t SharedMemorySize;

    //if (sizeof(float2) * THREADS_PER_BLOCK + SHARED_KERNEL<= props.sharedMemPerBlock)
		//SharedMemorySize = sizeof(float2) * THREADS_PER_BLOCK;
	//	SharedMemorySize = 1024;
	//else
	//{
	//	printf("sem shared: %d",(int)(sizeof(float2) * N + SHARED_KERNEL));
		//SharedMemorySize = props.sharedMemPerBlock - SHARED_KERNEL;
	//}

		//printf("\nTenho SharedMem : %d\nUsarei SharedMem: %d",(int) props.sharedMemPerBlock,(int) SharedMemorySize);
		//getchar();

	//if ( cudaSuccess != cudaGetLastError() )
    //{
	//	printf("\nErro antes do kernel: %s",cudaGetErrorString(cudaGetLastError()));
	//	getchar();
	//}

/*    EnergyCalculation<<<NUMBER_OF_BLOCKS, THREADS_PER_BLOCK, SharedMemorySize>>>( d_ParticleArrayPosition,
                                                                                  d_ParticleArrayVelocity,
                                                                                  N,
                                                                                  epson,
                                                                                  dt,
                                                                                  d_KineticEnergy,
                                                                                  d_PotentialEnergy,
                                                                                  trace,
                                                                                  d_traceKE,
                                                                                  d_tracePE );
   if ( cudaSuccess != cudaGetLastError() )
    {
		printf("%s line %d: CUDA Error: %s\n", __FILE__, __LINE__, cudaGetErrorString(cudaGetLastError()));
		getchar();
	    //printf( "Error kernel EnergyCalculation!\n" );
	}
*/

    // copy calculated energy values from device memory back to host memory.
   // CUDA_SAFE_CALL( cudaMemcpy(KineticEnergy,         d_KineticEnergy,         sizeof(REAL),    cudaMemcpyDeviceToHost) );
   // CUDA_SAFE_CALL( cudaMemcpy(PotentialEnergy,       d_PotentialEnergy,       sizeof(REAL),    cudaMemcpyDeviceToHost) );

    //teste
    //printf("\nParticleArrayPosition: x:%f y:%f",ParticleArrayPosition[0].x,ParticleArrayPosition[0].y);
    //CUDA_SAFE_CALL( cudaMemcpy(ParticleArrayPosition, d_ParticleArrayPosition,   sizeof(REAL),    cudaMemcpyDeviceToHost) );
    //printf("\nParticleArrayPosition: x:%f y:%f",ParticleArrayPosition[0].x,ParticleArrayPosition[0].y);

    //CUDA_SAFE_CALL( cudaMemcpy(traceKE,       d_traceKE,       sizeof(REAL),    cudaMemcpyDeviceToHost) );
    //CUDA_SAFE_CALL( cudaMemcpy(tracePE,       d_tracePE,       sizeof(REAL),    cudaMemcpyDeviceToHost) );

    REAL InitialKE = *KineticEnergy;
    REAL InitialPE = *PotentialEnergy;
    REAL InitialTotalEnergy = InitialKE + InitialPE;
    traceParticle(trace,t,traceKE,tracePE,"w");
    EnergyLog("./data/energy.log", 0, InitialTotalEnergy , InitialKE, InitialPE, 0.0, "w");


    while (t<=timeSteps)
    {

        // run the Kernel
        //UpdateParticlesAcceleration<<<1, N>>>( d_ParticleArrayPosition, d_ParticleArrayVelocity, N, 4.0/395.0, dt );
        //UpdateParticlesAcceleration2<<<NUMBER_OF_BLOCKS, THREADS_PER_BLOCK>>>( d_ParticleArrayPosition, d_ParticleArrayVelocity, N, epson, dt );

/*        UpdateParticlesAccelerationAndEnergyCalculation<<<NUMBER_OF_BLOCKS, THREADS_PER_BLOCK,SharedMemorySize>>>( d_ParticleArrayPosition,
                                                                                                                   d_ParticleArrayVelocity,
                                                                                                                   N,
                                                                                                                   epson,
                                                                                                                   dt,
                                                                                                                   d_flag);


__global__ void MainKernel(REAL_CUDA* d_ParticleArrayPosition,
                           REAL_CUDA* d_ParticleArrayVelocity,
                           unsigned int N,
                           REAL SofteningFactor,
                           REAL TimeStep,
                           REAL	*d_flag, // para que serve isso, Lucas??? Vou manter por enquanto.
                           int tileWidth)


*/                                                                                                                    

		// run the main kernel
		MainKernel<<<NUMBER_OF_BLOCKS, THREADS_PER_BLOCK,SharedMemorySize>>>(d_ParticleArrayPosition, d_ParticleArrayVelocity, N, epson, dt, d_flag, tileWidth);
       
       
       cudaError err = cudaGetLastError();
        if ( cudaSuccess != err )
        {
			printf("ERRO NO UPDATE... %s line %d: CUDA Error: %s\n", __FILE__, __LINE__, cudaGetErrorString(err));
			getchar();
			//printf( "Error kernel UpdatePaticulesAndEnergyCalculation!\n" );
		}


        if (t <= 0.0002)
        {
			CUDA_SAFE_CALL(cudaMemcpy(&flag,d_flag, sizeof(REAL), cudaMemcpyDeviceToHost));
			printf("\n-------- FLAG = %f",flag);
		}


        t += dt;

        if(t>tout)
        {
                // copy the updated data arrays from device global memory back to host memory.
                CUDA_SAFE_CALL( cudaMemcpy(ParticleArrayPosition, d_ParticleArrayPosition, sizeof(REAL_CUDA)*N, cudaMemcpyDeviceToHost) );
                CUDA_SAFE_CALL( cudaMemcpy(ParticleArrayVelocity, d_ParticleArrayVelocity, sizeof(REAL_CUDA)*N, cudaMemcpyDeviceToHost) );
               // CUDA_SAFE_CALL( cudaMemcpy(KineticEnergy,         d_KineticEnergy,         sizeof(REAL),    cudaMemcpyDeviceToHost) );
               // CUDA_SAFE_CALL( cudaMemcpy(PotentialEnergy,       d_PotentialEnergy,       sizeof(REAL),    cudaMemcpyDeviceToHost) );

               // CUDA_SAFE_CALL( cudaMemcpy(traceKE,       d_traceKE,       sizeof(REAL),    cudaMemcpyDeviceToHost) );
               // CUDA_SAFE_CALL( cudaMemcpy(tracePE,       d_tracePE,       sizeof(REAL),    cudaMemcpyDeviceToHost) );

                tout = t + dtout;
                sprintf(ResultFileName, "./data/%03d.dat", i++);

                printf("\n::: gerando arquivo de dados %s...\n", ResultFileName);
                GenerateResultFile(ResultFileName, ParticleArrayPosition, ParticleArrayVelocity, N);

                // Energy Logging
                REAL EK = *KineticEnergy;
                REAL EP = *PotentialEnergy;
                REAL EnergiaTotal  = EK + EP;
                REAL EnergyError = ( EnergiaTotal - InitialTotalEnergy) / fabs(InitialTotalEnergy);
                traceParticle(trace,t,traceKE,tracePE,"a");
                EnergyLog("./data/energy.log", t, EnergiaTotal , EK, EP, EnergyError, "a");

                somErr += EnergyError;
                contErr++;
                //EnergyLog("./data/energy.log", EnergiaTotal , 0.0, 0.0, "a");
        }

    }

    FILE* file_p;
    file_p = fopen("./data/errEnergia.txt","w");
    fprintf(file_p,"\n=======MEDIA ERRO ENERGIA=========\n");
    fclose (file_p);
    file_p = fopen("./data/errEnergia.txt","a");
    fprintf(file_p,"%e",somErr/contErr);
    printf("%e",somErr/contErr);
    fclose (file_p);


    // copy the updated data arrays from device global memory back to host memory.
    CUDA_SAFE_CALL( cudaMemcpy(ParticleArrayPosition, d_ParticleArrayPosition, sizeof(REAL_CUDA)*N, cudaMemcpyDeviceToHost) );
    CUDA_SAFE_CALL( cudaMemcpy(ParticleArrayVelocity, d_ParticleArrayVelocity, sizeof(REAL_CUDA)*N, cudaMemcpyDeviceToHost) );

    cudaFree(d_ParticleArrayPosition);
    cudaFree(d_ParticleArrayVelocity);
    free(ParticleArrayPosition);
    free(ParticleArrayVelocity);



    printf("\nValor de N: %d\n", N);


    printf("\n\nYou ran me by typing: %s\n", argv[0]);
    printf("\nPlaca: %s\n", props.name);
    return (0);
}
